Relativistic stars in f(R) and scalar-tensor theories 



Eugcny Babichev a,& , David Langlois a,c 
a ARC, UMR 7164 (G'NRS-Universite Pans 7), 10 rue Alice Domon et Leonie Duquet,75205 Pans Cedex 13, France 
b Institute for Nuclear Research of the Russian Academy of Sciences, 
60th October Anniversary Prospect, 7a, 117312 Moscow, Russia and 
c IAP (Institut dAstrophysique de Paris), 98bis Boulevard Arago, 75014 Paris, France 

(Dated: August 5, 2010) 

We study relativistic stars in the context of scalar tensor theories of gravity that try to account 
for the observed cosmic acceleration and satisfy the local gravity constraints via the chameleon 
mechanism. More specifically, we consider two types of models: scalar tensor theories with an inverse 
power law potential and f(R) theories. Using a relaxation algorithm, we construct numerically static 
relativistic stars, both for constant energy density configurations and for a polytropic equation of 
state. We can reach a gravitational potential up to $ ~ 0.3 at the surface of the star, even in f(R) 
theories with an "unprotected" curvature singularity. However, we find static configurations only if 
f ^ , the pressure does not exceed one third of the energy density, except possibly in a limited region of 

■ the star (otherwise, one expects tachyonic instabilities to develop). This constraint is satisfied by 

realistic equations of state for neutron stars. 

bjo; 

^ . PACS numbers: 04.50.Kd, 04.40.Dg, 95.36,+x 
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O" 1 I. INTRODUCTION 

U . 

One of the most challenging tasks for cosmology and fundamental physics today is to try to understand the apparent 
acceleration of the Universe. Beyond the minimal assumption of a pure cosmological constant, two main approaches 
£T) . have been explored. The first one consists in assuming some unknown form of matter, called dark energy, characterized 
^ ' by an equation of state P ~ —p. The second approach is more radical, as it tries to explain the present observations 
r**"* , as the manifestation of a modified theory of gravity, which mimicks general relativity on solar system scales, but 
ON ■ significantly deviates from it on cosmological scales. 

It turns out that it is rather difficult to construct a theory of gravity that, while being internally consistent, 
. ■ can account for the cosmological observations and be compatible with the present gravity constraints deduced from 
laboratory experiments and from solar and astrophysical systems. A class which has attracted a lot of attention is 
the so-called f(R) gravity theories where the standard Einstein- Hilbcrt gravitational Lagrangian, proportional to the 
scalar curvature R, is replaced by a function of R while the matter part of the Lagrangian is left unchanged (see e.g. 
[l[ for a recent review). After several detours, it has been realized that viable f(R) theories must satisfy stringent 
conditions in order to avoid instabilities and to satisfy the present laboratory and solar system constraints, and a few 
models have been carefully constructed to meet these requirements @-13|- 

To explore the full viability of these theories, it is important to go one step further by studying their behaviour 
in the strong gravity regime, such as reigns in the core of the most relativistic stars, namely neutron stars. In this 
context, it has been claimed in that very relativistic stars do not exist in the models because of the presence 
of an easily accessible singularity. In a recent work Q, we have shown that this claim does not hold by constructing 
numerically static relativistic stars. Note that relativistic stars have also been studied in but within different f(R) 
models. 

f(R) models can also be seen as a subclass of scalar-tensor theories. In particular, viable f(R) models rely on the 
so-called chameleon mechanism. For this reason, it is interesting to extend the study of relativistic stars to chameleon 
models. In this work, we show that the behaviour of chameleon and f(R) models is quite similar. In particular, 
in both cases, the scalar field in the innermost part of the star sticks to the minimum of its effective potential (if 
this minimum exists). If the equation of state is such that p — 3P < 0, which occurs in the central part of highly 
relativistic constant energy density stars, then there is no minimum. As we show here, it is nevertheless possible to 
construct numerically static stars, up to some critical value of the central energy density. We believe that this is due 
to tachyonic instabilities, associated with a negative effective squared mass, which develop and prevent the existence 
of a static star configuration. However, this problem does not apply to realistic neutron stars: although there is a 
large uncertainty on the equation of state deep inside a neutron star, the equations of state that have been proposed 
in the literature verify p — 3P > throughout the star. To construct a simple approximation of a realistic neutron 
star, we have used a polytropic equation of state. 
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Although, according to our previous work [7| and the present one, the singularity of the models (2j-|4| docs not 
seem so far to be an obstacle for relativistic stars, it appears to be problematic for cosmology p|. This has motivated 
the construction of regularized versions of these models by adding for instance an extra R 2 term 0, A more 

sophisticated model was also proposed recently in [l3[. We will consider these "cured" f(R) theories and compare 
their behaviour in relativistic stars with their "singular" counterparts. 

This paper is organized as follows. In the next section, we derive the main equations governing static and spherically 
symmetric configurations in scalar-tensor theories. Section III is devoted to relativistic stars in chameleon models. 
We then consider, in Section IV, the f(R) models, including the "cured" models that have recently been introduced 
to solve the cosmological singularity problem of some of these models. We finally conclude in Section V. 



II. STATIC AND SPHERICALLY SYMMETRIC EQUATIONS 



In the present work, we consider models characterized by an action of the form 



S 



d 4 Xy/—g 



(1) 



with Mp = 1/(8ttG), and where the matter (i.e. the fields <t m ) is minimally coupled to the metric 

g^ u = n 2 (cf))g^. 
For simplicity, we consider only couplings of the form 



£7 = exp Q 



dp) ' 



(2) 



where Q is a constant, although more general functions can have interesting effects (such as the spontaneous scalar- 
ization dicovered in 14J for coupling functions of the form lnfi = (3(j) 2 /2 with (3 < 0). 

Viable scalar-tensor theories are severely restricted by the present constraints on gravity, coming from laboratory 
experiments, solar system and binary pulsar tests. For example, deviations from General Relativity (GR) can be 
parametrized by the post-Netwonian parameters — 1) and (7 — 1). For the coupling ([2]), one finds (see e.g. [Till 1 ) 



/3-l = 0, 



7 - 1 = -4 



Q 2 



1 + 2Q 2 



The present constraint, inferred from solar system tests, is |16| 

I 7 - 1 |< 2 x 1(T 5 . 

The binary pulsars also give a similar constraint, but which is so far weaker than the solar system constraint in our 
case 2 . The above constraint suggests that scalar tensor theories of gravity can be viable only for very small couplings 
between ordinary matter and the scalar field. 

However, this constraint can be evaded if the scalar field is endowed with a potential such that its effective mass 
becomes large in the presence of matter. This is the so-called chameleon effect [ID, 111. As a consequence of this 
effect, the "bare" coupling constant Q is replaced by an effective coupling Q c s which can be strongly suppressed, 
i.e. such that |Q c ff| "C \Q\- The purpose of this paper is to study relativistic stars in this type of gravity. This also 
applies to f(R) theories, which can be recast in the form of scalar tensor theories with Q = — l/v6, leading to the 
unacceptable value 7=1/2 1201. Therefore, these theories can be considered as viable if the corresponding potential 
allows for a chameleon effect I2l|-|23ll. 



1 Note that our <j> is related to the scalar field tp defined in |15| by <j> = \/2Mpip and therefore their a corresponds to a = \/2Q. 

2 This can be seen on the figure 11 of the very recent lecture notes by G. Esposito-Farese |17| . For models with a coupling of the form f2l , 
the parameter space is restricted to the vertical axis /9o = where, so far, the solar system constraints turn out to be more stringent. 
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A. Links between the Jordan and Einstein frames 



The action ([T]) is defined in the Einstein frame, where gravity is described by the usual Einstein-Hilbert term for 
the metric g^ v . This metric differs from the metric (L„, defined in the Jordan frame, which is directly felt by the 
matter. As a consequence, one must be careful to specify the metric with respect to which various quantities are 
defined. For example, the energy-momentum tensor, defined in the Jordan frame and denoted T M „, is related to the 
energy- momentum tensor defined in the Einstein frame, T„„, by 

f£ = fr 4 T^, 

so that, for a perfect fluid, we have 

p = n 4 p, P = Q 4 P. 

The equation of state will always be specified in the Jordan frame. When \Q<fi/Mp\ <§; 1, as will be the case for our 
relativistic stars, the numerical values for the energy density and the pressure are essentially the same in the two 
frames. 



B. Equations of motion 

In Einstein's frame, the equations of motion are given by 



Guv = Ml 



v ff v c 



dV w 



rp(m) 



(3) 
(4) 



where G av = — (l/2)Rg„ I/ is the Einstein tensor, Tjffl is the energy-momentum tensor for the fluid matter of 
the star (i.e. does not include the scalar field) and 

T (m) = g^T^ = -p + 3P 

is its trace. 

We now consider a static and spherically symmetric geometry, with metric 

ds 2 = -e v di 2 + e x dr 2 + r 2 (d0 2 + sin 2 6 d(j) 2 ) . 
Introducing the radial function m(r) so that 

e~ A = 1 - 2m/r, 

the time and radial components of Einstein's equations ([3]) yield, respectively, 



2M 2 P 



o 4 p + - e -V 2 + n0) 



v' = e x 



M 2 P 



(5a) 
(5b) 



where one recognizes, inside the brackets of the first equation, the total energy density, which includes the fluid energy 
density (defined in the Einstein frame) as well as the gradient and potential energies of the scalar field. 

Instead of the angular component of Einstein's equations, it is convenient to use the conservation of the fluid 
energy-momentum, which has the standard form, 

in the Jordan frame, where the matter is minimally coupled. This gives 



(G) 
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The last equation is provided by the Klein-Gordon equation, Eq. for the scalar field, 

0"+0 + ^'-A')) 4>' = e x 

Finally, an equation of state, 

P = P{p), (8) 

closes the system of equations (l5al) . (|5b|) . (JB]) and ([7]). 



dV 



fl 3 n'(p - 3P) 



(7) 



C. Constant energy density stars 

In general relativity, one can solve analytically the profile of a relativistic star by assuming that the energy density 
is constant, 

P = Po ■ 

Using (|5a)) in the GR limit, which corresponds to <j> = and V{4>) = 0, this implies 



6Mf, 



Substituting in the Tolman-Oppcnheimer-Volkov equation 

p + P 



P' 



1M 2 p 



r 2 (l — 1m It) 

which follows from © and (|5b[) in the GR limit, one finds, after integration, that the pressure profile is given by 

2GM, \ 1/2 2GM,r 2 \ 1/2 



P M = \ 1/2 7 ' \l/2 ( 9 ) 



(l_2£^i) _ 3 (l-^) 



where r* is the radius of the star and M* = 47rpo^*/3 its mass. 

We can use the above analytical results to compute the trace of the energy momentum-tensor, which is the quantity 
that couples directly to the scalar field in its equation of motion. The relation © implies 



/ 9 \ 1/2 / \ 1/2 

4 (l- ^GM^ _ 6 (l- 2£M.) 

P ~ ^ = P °~ ^ T ;£ ,1/2 • 



(1-2^) _ 3 (i-2£M,y 

Since this is a function that grows with increasing r, its minimal value is at the center of the star and is given by 

1 -3(1 - 2$*) 1/2 

where $^ is the gravitational potential at the star surface. Therefore, whenever the compactness of the star is higher 
than the critical value 

*, = A ~ 0.28, (10) 

the quantity p — 3P becomes negative in the most central layers of the star. As we will see later, the sign of 
p — 3P is crucial to understand the behaviour of the scalar field inside the star and can have drastic implications 
concerning the stability of the scalar field. Note that this critical value was also pointed out in [24 |. Note also that 
although our derivation above applies only to constant energy density stars in pure general relativity, it remains a 
good approximation for the stars we consider because the backreaction of the scalar field on the star profile is very 
small. 
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D. More realistic equation of state 

Although analytically simple, a constant energy density star is not a very realistic substitute for a real neutron star. 
Since the equation of state deep inside a neutron star is still unknown, we use in this paper a polytropic equation of 
state, which remains simple and is believed to represent a reasonable approximation for a real neutron star. Following 
[25| we use, 

n + K—\, P(n) = Km B — , (11) 
n J n 

with raj = 1.66 x 1CP 27 kg, no = 0.1 fm~ 3 and K = 0.1. We take for the central particle number density the value 




FIG. 1: Energy density p (solid blue line), pressure P (dashed purple line) and the combination p — 3P (black dotted line), in 
units of the central density p c , as functions of the radial coordinate r (in units of Mp/5^~ 1//2 ). 

h c = 0.3 fm -3 , which corresponds to a central energy density p c = 6.47 x 10 14 g/cm 3 . In general relativity, the 
numerical integration of the Tolman-Oppenheimer-Volkov equations leads to a star, which is stable, with a global 
mass M = 3.10Af© and a gravitational potential |$*| ~ 0.21. In both the chameleon and f(R) models, we will use 
the same equation of state with the same central density and will find that the star configuration is almost identical 
to the general relativistic one. 

In Fig. [I] we plot the radial profile of the energy density and of the pressure. As can be seen in the same figure, 
the quantity p — 3P remains positive throughout the star. The same holds for other realistic equations of state that 
have been proposed in the literature (see e.g. (26j and references therein for a discussion on the equations of state of 
neutron stars) . Interestingly, p — 3P is not a monotonous function of the radius in the example we have chosen. We 
will see that this has a visible influence on the profile of the scalar field inside the star. 

E. Effective potential 

Assuming some configuration for the relativistic star and ignoring any backreaction of the scalar field on this 
configuration, it is useful to interpret the right hand side of the Klein-Gordon equation ([7]), up to the metric component 
e A , as the derivative of an effective potential defined by 



V cS = V+^Q 4 (p-3P). 



(12) 
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Consequently, a local extremum of the potential is characterized by the condition 

dV eS dV , 



^-e iQ *l Mp {~p-?,P) = 0. (13) 
Mp 



A solution of the above equation, if it exists, corresponds to a minimum of the potential if the effective squared mass 

M P 2 



= + 4-^e 4c ^ (p 3P) 



is positive. It is thus manifest that the quantity p—3P plays a crucial role in the overall sign of m 2 s , with consequences 
for the stability of the star, which we discuss just below. 



F. Stability 

Although the study of the stability of our star configurations would require a detailed analysis, which we leave 
for future work, one can try to discuss this issue with some qualitative arguments. For simplicity, we ignore the 
perturbations of the star matter and of the scalar field and consider only the perturbations of the scalar field about 
the background configuration. Decomposing these perturbations into spherical harmonics, 

6(f>(t, r, 6, ft = ^ 6<h m (t, r)Y lm (9, 0), 

one finds that each mode verifies the equation 

5<f) - e 

where m 2 s has been defined in p4[) . Using a Fourier decomposition in time, and denoting the frequency u>, the above 
equation of motion leads to a dispersion relation of the form 

2 ; 2 , 2 

where A: -1 is the typical lengthscale of variation of the scalar field. One can write k = ar" 1 , where a is typically 
bigger than 1 but not very much, in the situations we consider. Let us also assume that the dominant contribution 
to m 2 s comes from the matter part. Concentrating on the radial modes (which are potentially more dangerous), we 
find that the configuration is stable if, 

2 a 2 Q 2 , 

w ^7f + 4 Mf^ (1 ~ 3u;)>0 ' 

where 

w = P I p. 

Substituting p* ~ M*/r^ and $* ~ M*/ \Mp 2 r*) leads to the condition 

(1 - 3u))$* > -1, 

where the right hand side is some constant, but typically of order —1. According to this relation, obtained in a very 
crude way, one expects instabilities to appear whenever the equation of state is such that 1 — 'Sw < and the star is 
very compact so that $* is a significant fraction of one. This is precisely the situation which one encounters for very 
massive stars with constant energy density, where w > 1/3 in the central region. 



1(1 + 1) 



'off 



0. 



G. Numerical procedure 



We now present how we proceed numerically As a fist step, it is convenient to rescale the various quantities involved 
and to work directly with dimensionless quantities. Let us thus introduce the following rescaled variables 

M 2 M 2 M 2 

r = r Q Z, <t> = M <t>, V = -^-V, p = ^p, P* = ^P, (14) 
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where r 1 , Mo, M\ are parameters of dimension of mass, which are so far arbitrary. Substituting the above rescalings 
(TH|) in the equations of motion (|5aj) , (|5b[) . © and (J7J, we obtain the following system of equations, 



2M' 



-e" V' 2 +V), 



2M 



4Qe « 



v , 



P> = --[p + P) [v' + 2Qe (f>' 



O = 0"+(| + i(^-A') 



(15a) 
(15b) 
(15c) 
(15d) 



where we have used the particular coupling (|2|) and introduced the notations, 



M = -^{l-e 



Mq 
M P ' 



Mi 
Mp 



The values of the parameters eo and €i have no physical relevance. They are chosen purely for convenience. In the 
following, we will take e\ = 1 in all cases, but our choice for eo will be different for chameleon models and for f(R) 
gravity. 

Meanwhile the equation of state (jSJ transforms into 

P = P(p). (16) 
For the polytropic equation of state Eq. ([TTj) . with h c = 0.3 fm -3 , we will always choose ro = Mpp c 1 ^ 2 , so that 

p c = 1 (polytropic star). 
The equation of state Eq. ([TTj) then becomes, in these rescaled units, 



P = P 



P 
0^39 



P = 



1 



1.56 



y/l + 1.56/5 - 1 



(17) 



with the rescaled central pressure P c ~ 0.23. Since the neutron star is a highly relativistic object, this choice entails 
that the radius of the star in rescaled units is of order one, = r*/ro ~ 1. 

Since we study static configurations, a natural way to proceed with numerics is to use a relaxation algorithm. In 
the following calculations the relaxation parameter will be taken to vary with the iteration number (typically between 
10~ 3 and 1), in order to ensure the convergence of the algorithm and to minimize the number of iterations. We 
have used a non-homogenous grid, usually between 1000 and 6000 points. The details of the grid depend on the 
model under consideration and its parameters. Normally the grid is chosen so that the points of the grid are highly 
concentrated in the particular region inside the star where the solution varies rapidly. Such a choice of inhomogeneous 
grid allows us to numerically resolve potentially problematic regions. 

Let us finally discuss the boundary conditions. The system (|15[) . which includes three first-order differential equa- 
tions and one second-order differential equation, requires five boundary conditions. Some of them are defined at the 
center of the star, at £ = 0, while the remaining ones are specified far from the star, at some point £ = £2- The 
boundary conditions for the metric components and the scalar field are 



M(o) = o, K6) = o, ci'(o) = o, 0(6) = k 



(18) 



The boundary conditions for M. and <f> at the center just follow from the requirement of regularity at £ = 0. The 
condition on v is in fact arbitrary, since only v' enters (|15|) but not v itself: our choice corresponds to defining the time 
coordinate as the proper time of a static observer at £ 2 - The "asymptotic" value corresponds to the minimum of 
the effective potential V e s far from the star 3 . Finally, the boundary condition for the pressure is 



P(0) = P c 



Some situations require a more precise boundary condition, in which case the exact value of <j> at £2 can be deduced from the asymptotic 
behaviour of the scalar field at large distances, which can be computed analytically. 
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depending on the problem: for a constant-density star it is convenient to choose P(&) = —poo so that the solution is 
asymptotically de Sitter, while for the polytropic equation state, the pressure at the origin is fixed 4 . 

III. CHAMELEON 
A. Description of the models 

In this section we consider chameleon scalar-tensor theories, introduced in [l8l fl9j. In these models, the effective 
mass of the scalar field becomes important inside matter so that the usual constraints on fifth-force interactions are 
evaded. To illustrate the behaviour of a chameleon field in a relativistic star, we will concentrate on a specific model, 
but our results should remain qualitatively valid for other models. 

We will consider the action ([T]) with the potential 

V&) = 

where p is a mass parameter. We will also assume that Q is positive. Moreover, as in the original proposal for the 
chameleon 0, we will use Q ~ 1 and p ~ 10~ 3 eV (this value corresponds to an upper bound for p in order to 
satisfy the equivalence principle constraints [3). 

If Qcp/AIp <C 1, which will be true for our configurations, the scalar field value corresponding to the minimum of 
the effective potential (fl~2j) is given by 



Mpfi 5 



Q(p-3P) 



Note that this is defined only if p — 3P > 0. Therefore, there is no local minimum at the center of a constant energy 
density star with > 5/18. 

For our numerical study, we have chosen eo such that the asymptotic value of the scalar field is normalized to a 
value of order unity. More precisely, we take 5 



eo = 



QMppoo 

where poo is the asymptotic energy density, so that the local minimum is given by 



(1 - 3w)p ' 

in rescaled units. Asymptotically, one thus finds </>oo = 1/^/1 — Sw^ = 1/2 (for Woo = — 1). 

As for the parameter ro, it is convenient to choose it of the order of the star radius. In the case of the polytropic 
equation of state, we take ro = Mpp c , so that p c = 1, as mentioned earlier. For the constant energy density stars, 
we will adopt the same prescription for a reference star characterized by $* = 0.165. We will then compare this star 
to other stars, with the same radius but different central densities, by keeping ro fixed but by allowing the rescaled 
central density p c to vary. 

The rescaled potential can be written in the form 

V = — , V = — poo . (19) 

(f) £o 

In rescaled units, taking ei = 1, the effective potential is thus 

Vcs = -^e 4 ^05-3P). (20) 
6 4 eo 



4 In this case we also obtain asymptotically de Sitter solution. 

5 Any other choice is acceptable, but not so well adapted for numerics. For example, the choice eo = 1 leads to values of the rescaled 
scalar field that are very small fractions of unity, which is not very convenient for the numerical investigation. 
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while the corresponding effective squared mass is 

- 2 2v o 



'eff 



>(p-3P). 



(21) 



Let us estimate the order of magnitude of the above parameters. A typical value for the density inside neutron 
stars is p c ~ 10 15 g/cm 3 , whereas the asymptotic density is p^ ~ 10 _24 g/cm 3 , corresponding to the galactic density. 
With Q ~ 1 and /i ~ 10~ 3 eV, this leads to tiny values for our dimensionlcss parameters: eq ~ 10 -19 and vo ~ 10~ 21 . 
These realistic parameters are too small for the numerical precision that can be reached in practice and we have used 
parameters that are much bigger. 




FIG. 2: The rescaled bare potentials (|19|) . shown by solid blue, the matter part of the effective potential, shown by dashed red, 
and the effective potential for the chameleon model (|20[) . shown by dotted black line. The parameters are chosen as follows, 
Q = 1, vo = 0.05, eo = 0.05, ei = 1 and p - 3P = 0.14. 

Asymptotically, the expression yields (|2"Tj) the tiny value 

m 2 ff:OO ~2i;„(l-3«; 0O ) 3 / 2 , 

where the second term is negligible for eo <C 1. By contrast, at the center of the star, the minimum is extremely small, 

^ min < c = \lQ(l-3w c y (22) 

while the effective mass becomes huge 



m eS,c - 2 \ 



IQ 3 (l-3w c ) 3 
v o4 



One can thus expect that the scalar field configuration will follow its local minimum (when it is defined) in the denser 
regions of the star. 



B. Numerical results 



Using the full system of equations combining Einstein's equations and the Klein-Gordon equation, we have computed 
numerically the profile of the chameleon field, of the matter and of the geometry for two types of star configurations. 
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Before discussing these two cases, let us comment about a subtlety concerning the asympotic behaviour of the 
geometry far from the star. In a realistic context, the star configuration should be matched asymptotically to a 
cosmological geometry, which is time dependent. Since we restrict our analysis to purely static configurations, we 
introduce instead some artificial matter which behaves like a cosmological constant far from the star. At large radius, 
the geometry thus approaches a Schwarzschild-de Sitter metric, which is static. We have checked that the details of the 
asymptotic matter do not modify the star configuration. Another possibility would be to introduce some non-trivial 
minimum in the chameleon potential. 



1. Constant energy density stars 

As a first (toy) example we considered a constant energy density star. In order to avoid unnecessary numerical 
difficulties connected with a sharp transition of the density from the star to the surrounding medium, we introduced 
the following smoothed profile for the density: 



HQ = f 



1 — tanh 



Poo, (23) 



where p c is the rescaled density at the center of the star, p^ is the density of the surrounding medium far from the 
star and a is a "smoothing" parameter (with values between 10 -12 and 10 -2 ). It is easy to see, that for small enough 
a the density is approximately equal to p c for £ < 1 and p ~ p^ for £ > 1 , so that £ ~ 1 is the radius of the star. 
The boundary condition for the pressure is, P{^2) = ~Poo, so that we ensure that the de Sitter asymptotic behaviour 
is recovered for our numerical solution. In this case the pressure is found as a function of £, and P inside the star 
depends on the boundary condition £2- This seemingly paradoxical situation is due to the artificial choice (|23|) . but, 
as we have checked, the effect is tiny because the ratio between the star energy density and the asymptotic one is 
huge. 

For a given coupling Q, the numerical system depends on three independent parameters: p c , vq and eo, while 
4>oo = 1/2 and the rescaled asymptotic density poo is determined by poa — vq£o/Q according to (|19[) . Since the 
backreaction of the scalar field on the star itself is very small, the compactness of the star depends only on the 
parameter p c . Consequently, for fixed parameters vq and eo, increasing the central density leads to stars that are 
more and more compact, and therefore more and more relativistic. In Fig. [31 we plot the profile of the scalar field 
for stars with different central densities, corresponding to a gravitational potential ranging from 0.00168 to 0.165. 
As the central density increases, one sees how the scalar field profile evolves from a smooth configuration to a thin 
shell configuration, where the scalar field abruptly jumps within a small range of radii. This is in agreement with the 
Newtonian result [lj] that Ar/r* « (0oo - c )/(6QM P $ >t ) « e /(l2Q$*). 

Essentially the thin-shell solution can be understood as follows. The solution tries to minimize the sum of the 
potential energy and of the gradient energy, with the constraints that 4>' = at r = and <f> fixed asymptotically. As 
the energy density in the star increases, to be away from the minimum of the effective potential becomes very costly 
energetically and, therefore, the field prefers to stay as much as possible close to its minimum. At some point close 
to the star radius, however, the field starts to move away from its effective minimum. In this region, the contribution 
from the potential on the right hand side of the Klein-Gordon equation becomes negligible with respect to the matter 
dependent term, and the field "evolves" very quickly, driven by this "matter force" term. In the relativistic regime, 
analytical expressions describing the profile of the scalar field in the thin shell regime have been obtained in [13] : they 
apply to constant energy density stars and in the linearized approximation for the gravitational potential. 

If one continues to increase the central energy density of the star, one ends up reaching the critical value for which 
p— 3P < in the innermost central region of the star. The profile of the scalar field in the central region of the star is 
then radically altered, as one can see on Fig. 01 because there is no local minimum for the scalar field in this region. 



2. Thin shell effect 

It is also instructive to see how much, depending on the density of the star, the geometry outside the star is close 
to that expected in general relativity. Numerically, one can estimate the post-Newtonian parameter 7 by comparing 
the coefficients of the terms in l/r in the metric components gtt and g rr . 

Since we work in the Einstein frame, let us first discuss the metric components g u and g rr . One subtlety here is that 
the metric is asymptotically Schwarzschild-de Sitter rather than simply Schwarzschild and the asymptotic behaviours 
of the metric components are therefore 

\ , 2GM A 2 

e" w e~ x fn 1 -r 2 . (24 

r 3 
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FIG. 3: Profiles of the scalar field <j> (in Planck units) as a function of radius of the star (in units of ro = Mppcpc 1 ^ 2 ) for constant 
density stars with the equation of state (|23|) . The parameters are chosen as follows, a = 0.01, Q = t\ = 1, eo = i>o = 0.01, and 
the rescaled densities of the star are (from light gray to black), p c = 0.01, 0.02, 0.05, 1. The values of the gravitational potential 
at the surface of the star are, respectively, 0.00168, 0.00334, 0.0083 and 0.165. The increasing rescaled density corresponds to 
an increasing physical density while the physical radius of the star is being fixed. 



To estimate M numerically, we have substracted from the metric components the cosmological constant term with 

where is the asymptotic minimum of the scalar field and Voo = V^ipoo). We have checked that the two metric 
components, in the Einstein frame, provide the same value M . This means, in particular, that the backreaction due 
to the scalar field gradient is negligible outside the star. 

In the Jordan frame, the metric components are g^ v = cxp(2Q(f>/Mp)g flv . The asymptotic behaviour of the scalar 
field is of the form 



GM 



effe 



where Q c s is the effective coupling of the scalar field to the star, which is strongly suppressed in the thin shell regime. 
Combining this expression with (|24|). it is easy to get 



9tt 



9rr 



2Q<t> 



1-2(1-2QQ, 



GM A 



oil' I 



1-2(1 + 2QQ cff 



r 

r 3 

GM A 



— r 
3 



where, for simplicity, we have neglected the exponential decay. The effective 7 parameter is therefore 

1 — 2QQ c ff 



7 



1 + 2QQ cff 



(25) 



The general relativistic value, 7 = 1, is thus recovered in the limit where the effective scalar charge Q e ff of the star is 
strongly suppressed, i.e. in the thin shell limit. This is illustrated in Fig. [5j where we have estimated numerically the 
effective coupling Q c g as well as the post-Newtonian parameter 7 for stars with increasing energy density. 
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FIG. 4: Profiles of the scalar field (f> (in Planck units) as a function of radius of the star (in units of ro = Mppcpc 1 ^ 2 ) for stars 
with p — 3P < in the central region. The rescaled densities of the star (from light gray to black): p c — 1, 1.7, 1.8, 1.9. The 
other parameters are the same as in Fig. The values of the gravitational potential at the surface of the star are, respectively, 
0.165, 0.280, 0.298 and 0.318. 

3. Polytropic relativistic stars 

Let us now consider relativistic stars with the more realistic equation of state (flT|) . Although this particular equation 
of state is appropriate for the bulk of the star, it must be slightly modified at the edge of the star for the following 
reasons. First of all, the equation of state ([XT)) does not describe the exterior of the star. However, in the case of 
the chameleon field it is crucial to have a non-zero background density, therefore for small densities Eq. ([T7]) must 
be changed to include the background (galactic or extragalactic) matter in the model. Moreover Eq. ([17]) leads to a 
problematic behavior in the limit p — > (i.e. at the edge of the star), where the equation of state becomes P oc p 2 . 
In this case, the relativistic Euler equation become ill-defined. One can also see this problem when describing the 
dynamics of a perfect fluid with P oc p 2 in terms of a scalar field: the emergent metric for perturbations becomes 
singular when p — > 0, as it is explained in detail in [28|. We choose the following modification of the equation of state 

o, 

where we have introduced the (rescaled) asymptotic energy density p^, which is an extremely small parameter with 
respect to the pressure and energy density in the star. It should be stressed for P > poo the modified equation of 
state ([26]) takes the original form of the equation of state (fl7|) . while for small P the modified equation of state allows 
for the cosmological term behavior 6 . Thus the modification of the equation of state only affects a tiny region at the 
edge of the star. We also tried other modifications for the equation of state, and we obtained the same results, up to 
small differences around the edge of the star, as expected. It is easy to see that for small p^ the "fluid" (|2l))) is in the 
regime of a cosmological constant when — P = p ~ poo. 



In fact, our numerical integration shows that the system settles down at the cosmological term, P = 



— p, at large £. 
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Having specified the equation of state for the star and the surrounding matter (|26|) , we also need to fix the boundary 
conditions. Apart from the four usual conditions (|18[) , the boundary condition on the pressure is specified at the center 
of the star: P(0) ~ 0.23 (which corresponds to a realistic value, in the rescaled units). It is interesting to note, that the 
system "chooses" by itself the cosmological term-like behavior asymptotically, P = —p, we do require any conditions 
on the density nor the pressure at infinity. This situation is in contrast to the constant density star configuration, 
where we imposed the "correct" asymptotic behavior at infinity specifying the boundary condition for P. 

The numerical solution for the scalar field profile is shown in Fig. |6]for the parameters Q = e\ = 1, cq = 0.01, 
t>o = 0.01, po = 10~ 4 . The minimum of the effective potential is indicated by a dashed line. In the central part of the 
star, the scalar field strictly follows this minimum. Note that, as the radius increases, the minimum first increases and 
then decreases. This non-monotonous evolution is the consequence of the non-monotonous evolution of the quantity 
p — 3P, which was pointed previously. Slightly before reaching the radius of the star (delimited by the quasi-vertical 
dashed line), the scalar field starts to move away from its minimum, with a steep gradient, and asymptotically evolves 
toward its asymptotic minimum far from the star. We have also studied numerically configurations with other sets of 
parameters, namely in the range eo ~ 0.1 — 10 -3 , vo ~ 0.1 — 10~ 3 , po ~ 10~ 2 — 10~ 4 , and found a similar behaviour, 
while the value of the scalar field at the center of the star agrees with the estimate (f2"2"| . 

As mentioned earlier, the profiles for the energy density and pressure of the star are almost unchanged with respect 
to those for the corresponding star in pure general relativity (i.e. without scalar field), presented in Fig. [T] This can 
be understood by noting that the contribution of the scalar field in the first three equations (fT5|) can be neglected, 
giving then the same equations as in general relativity. The change in the total neutron star mass is less than one 
percent. 



IV. f(R) GRAVITY 

We now turn to f(R) gravity, restricting ourselves to the conventional metric formulation. Some of our results on 
relativistic stars have already appeared in a separate publication Q. Here, we present more details and stress the 
strong analogy of these models with the chameleon models. 




FIG. 5: Evolution of the numerical estimates for Q c s/Q (red triangles) and for the post-newtonian parameter 7 (blue dots) 
as we increase the energy density of the star (constant energy density star). The "bare" coupling is here Q — 1/2, which 
corresponds to 7 = 1/3 in the absence of screening. The other parameters are eo = 1(T 2 , v = 1CT 8 and £* = 1. 
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FIG. 6: Profile for the scalar field <j> (in Planck units), shown by solid (blue) line, as a function of the radius (in units of 
ro = Mpfic 1 ^ 2 ), for the equation of state (|26jl with eo = 0.01, vo = 0.01 and poo = 10 -4 . The value of the gravitational 
potential, .M (£)/£, at the surface of the star is 0.211. The value m i n for the minimum of the effective potential is plotted by 
dashed (gray) line. 



A. The models 

The f(R) models are usually described by an action of the form 

M 2 f 

S =^fj d*x^f(R) + Smi^g^}, (27) 

where the matter is minimally coupled to the metric g^vi with corresponding Ricci tensor R^ v and scalar curvature 
R. It is convenient to reexpress / in the form 

f(R) = RqT{x), x=^- 

where J 7 is a dimensionless function and Rq is some parameter. 

In contrast with many papers on f(R) theories, in particular Q and 0, we choose here to reexpress the model 
as a scalar- tensor theory in the so-called Einstein frame. The two formulations are of course equivalent (at least at 
the classical level), but the Einstein frame is useful to compare directly the behaviour of f(R) theories with that of 
chameleon models, discussed in the previous section. 

By introducing the scalar field 



and the metric 



^-M P \nf tR = sp-M P \nP{x) 



g^ = O 2 g^, Q 2 = / b = expU/ \ <f>/M P ] 



(28) 



the action ([2T]) can be reexpressed in the form ([T]) with the coupling 



Q = -1/VQ, 
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and the potential 



V 



M 2 P 



Rf,H-f 



(29) 



which can be expressed in terms of <p by inverting the definition (|28j) of <j) as a function of R. 

Note that, in contrast to the chameleon case discussed in the previous section, the coupling Q is here negative. 
We could have chosen a positive value by changing the sign of </>, but we have kept the usual convention. When 
comparing the results of this section with the previous one, it is useful to keep in mind that the signs of Q and <f> can 
be simultaneously changed. 




FIG. 7: Potential V (in units of MpRo) as a function of <f) (in Planck units) for n = 1 and = 3.6. The lower black dot 
corresponds to the de-Sitter attractor while the upper-right dot shows the curvature singularity. 



The function f(R) is severely constrained by observations, 
the specific class of models introduced by Starobinsky 



For the present discussion, we focus our attention on 



f(R) = Ro 



A l 



R 
Rq 



(30) 



The asymptotic de Sitter solution R^ = x^Rq corresponds to a minimum of V and depends on the parameter A. 
For practical purposes, it is simpler to choose x^ as the parameter and then express A as a function Xoo, with the 
restriction that Xoo must be chosen such that it corresponds to a stable minimum of the potential. Details and explicit 
expressions can be found in the appendix. The potential for the scalar field in the model with n = 1 and = 3.6 is 
plotted in Fig. 

B. High curvature regime 

Inside the star, the curvature is much higher than the cosmological curvature Rq, i.e. i» 1. It is instructive to 
consider the asymptotic expansion of the potential and its derivatives in order to understand intuitively the behaviour 
of the scalar field, noting that the infinite curvature limit corresponds to — > 0~. 

Using the asymptotic expansion 



T(x) 



X + \x 



-2n 



(X » 1), 



(31) 
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one finds 

V(x) = ~ - [1 - (1 + 2n)x~ 2n ] , 

which shows that the amplitude of the potential remains finite when x 3> 1, i.e. when <fi approaches 0. This is 
quite different from the chameleon model discussed in the previous section and this can cause some worry that the 
singularity can be easily accessible for high enough curvatures, as discussed in Q. 

However, the derivative of the potential still goes to infinity in this limit and this property prevents the scalar field 
to reach the singularity. Indeed, the derivative of the effective potential is given by (see appendix) 

— 5* ~ J-MpR - - -= (p - 3P , 

d4> V3 u 2 y/&M P K ' 



which implies the existence of a minimum 



f. < </>min = -V&Xr 



nx^- x . (32) 



M P 2 R 



at least if the matter term p — 3P is positive. It can be checked that this is indeed a minimum by computing the 
second derivative of the effective potential. This leads to a positive effective square mass (see appendix), 



2n+2 
/ > min 



6An(2n+ 1) 

which becomes very large at high curvature. 

As the energy density in the star increases toward the center, one expects that the scalar field will be closer and 
closer to the singularity (without reaching it however), simply because the minimum of the effective potential is closer 
and closer to the singularity. This is confirmed by our numerical calculations, which we dicuss below. 



C. Numerical results 



We have integrated numerically the system of equations ([T5]) . together with (fTTj) . for Starobinsky's model (f3"U)) with 
n = 1. Our parameters are eo = e% = 1 and the rescaled potential is 

V = v V, Vo = rtR = ^t^, 

Pc 

where the explicit expression for V in the case n = 1 is given by (| A6|) in the appendix. 

Here, the parameter w corresponds to the ratio between the energy density at infinity (i.e. the cosmological energy 
density) and the energy density at the center of the star. Realistic values of this parameter are thus extremely 
small, typically vo ~ 10~ 40 , and are numerically challenging to reach, because the scalar field value at the center is 
proportional to Vg, as can be seen by comparing the definition of vq with the minimum defined in (|32p in the case 
n=l. We have performed our numerical calculations in the range vq ~ 10 _1 — 5 x 10 -5 and checked that the solution 
at the center behaves as expected from our analytical analysis, Eq. (f3"2"]> . We are thus confident that this should hold 
for smaller values of vq. As in the case of the chameleon model, the profiles for the density and pressure only slightly 
differ from the GR results. 

The profile of the scalar field inside and outside the star is plotted in Fig. [S] As is clear from the figure, the scalar 
field tends to interpolate between an extremely high density regime, inside the star, and a very low density regime, 
outside the star. This behaviour is quite analogous to that of the chameleon model discussed in the previous section 
and the solution that we obtain is analogous to the thin-shell solution. What differs from the usual chameleon models 
is the presence of the singularity. In the very high density regime, the scalar field is very close to the singularity 
because the minimum of the effective potential is itself very close to the singularity. To give an analogy, it is like the 
scalar field is following a track very near a precipice, but because the effective mass is very high in this regime, the 
scalar field is securely attached to the track and does not fall into the nearby precipice. 

These numerical results, previously presented in Q, contradict the results of @, where it was claimed that stars 
with a gravitational potential larger than a critical value $ max m 0.1 could not be constructed. The reason advocated 
in Q was the presence of the singularity, as stressed previously in However, as we showed in Q, the analytical 
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FIG. 8: Profile of the scalar field 4> ( m Planck units), shown by solid (blue) line, as a function of the radius (in units of 
M P p~ 1/2 ), for the model (O with n = 1, a^oo = 3.6 (see appendix for the definition of !„) and no = 10 4 . The value </> m i n for 
the minimum of the effective potential is plotted by dashed (gray) line. 



arguments in [5j and [6| supporting this interpretation were not valid 7 . We also indicated that the exploration of 
the scalar field profile was numerically challenging for small values of vq, because the distance between the scalar 
field value and the singularity is proportional to Vq. Since the value used in Q was vq ~ 10~ 6 . one could strongly 
suspect that the results obtained in [6[ were probably due to a numerical instability. This explanation was confirmed 
subsequently in p9^ . which tried to reproduce the numerical analysis as @, using the same technique, namely the 
shooting method. Indeed, they were able to construct highly rclativistic stars, but with a parameter vq ~ 10~ 2 . 

Although we show only our numerical results in the case of the polytropic equation of state, we have also performed 
similar numerical integrations for constant energy density stars. We have obtained quite similar results as in the 
chameleon model. For instance, we have been able to obtain numerical solutions up to the value = 0.33 for the 
gravitational potential 8 (with the parameters ej = e a = 1, v a = 0.01). 



D. "Cured" f(R) gravity 

Although the existence of neutron stars is, so far 9 , compatible with the f(R) models we have just investigated, the 
situation is not the same in cosmology. It has indeed been noticed in [3, [30, HH that the behaviour of these theories 
is pathological when one considers their past cosmological evolution, with the generic presence of a singularity in the 
not so distant past. As a way to "cure" this pathology at high curvature, it has been suggested to add R 2 term 10 . 
This extra term radically changes the behaviour of the potential at very high curvature, as illustrated in Fig. [JJ] for 



7 Note also that, obviously, the problem encountered in [f| cannot be explained by the sign of the trace of the energy-momentum tensor 
since their gravitational potential is well below the critical value 1 1101 . 

8 Note that a similar value for the maximum gravitational potential, <&* = 0.345, was obtained in [29| . 

9 It would be also worth considering a supernova explosion and the birth of a neutron star in the context of f(R) theories. The present 
work deals only with static solutions and not how they could be reached in an evolutionary process. 

Another way to avoid a pathological behavior would be to construct f(R) model separating the singularity with infinite energy barrier, 
as it was done in |32(I . It was argued, however, that this particular model does not satisfy local gravity constraints and SDSS data l33t . 
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the model 



f(R) = Ro 



a(i- (i + x 2 y 



a x 



R_ 

Ro 



(33) 



One sees that the corresponding potential has no singularity at <fi = and the scalar field can now take positive values, 
corresponding to high curvatures. 

The extra R 2 term has the additional advantage that it can drive inflation in the early Universe (see e.g. [HI, HH), as 
in the pioneering Starobinsky model (36[. Recently, however, it was argued in [l3[ that one cannot relate smoothly the 
inflation epoch with the present acceleration epoch with the above model (|33[) . As an example of a viable transition 
between two accelerating solutions, the model 



f(R) = (1 - c)R + ceki 



cosh(_R/e- b) 
coshfr 



R 2 

6M 2 ' 



R(} 



b + ln(2cosh&) 



(34) 



with < c < 1/2, was proposed in [l3j], where the new mass scale M must be of the order 10 12 GcV to reproduce the 
standard predictions of inflation. In the high curvature regime (i? 3> Ro), one finds 



f(R) ^ R-cRo + cee 2b e- 2 "/ e 



R 2 
6M 2 



The interior of a neutron star corresponds to the high curvature regime, but the curvature is many orders of magnitude 
smaller than M 2 . In this intermediate regime Rq <C R <C M 2 , the scalar field behaves like 



-2ce 2fc e- 2 ^ e 



R 

3M 2 



(35) 



The two terms on the right hand side are extremely small. However, the first term decays exponentially as the 
curvature increases and the second term will typically dominate at the center of a neutron star. 

We plot on Fig.[in]the profile of the scalar field in this model with small values for vq and a = Rq/(6M 2 ), although 
not small enough to be realistic. Once again, we observe that the scalar field follows the minimum of its effective 




FIG. 9: Potential V (in units of Mj,R ) as a function of 4> (in Planck units) for n = 1 and x x — 3.6. The blue curve corresponds 
to the original Starobinsky cosmological model with the curvature singularity at <j> = 0. The red curve represents the "cured" 
Starobinsky model with a = 1CP 3 . 
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FIG. 10: Profile of the scalar field in the model ([34]) with the parameters cr = 2 x 10 6 , c = 0.4, 6 = 4, v = 10 3 . 

potential in the core of the star. Near the boundary of the star, it deviates from the local minimum and converges 
towards the asymptotic minimum outside the star. The scalar field is negative deep inside the star because the second 
term on the right hand side of ()35j) dominates, whereas it is positive near the boundary and outside the star. 

V. DISCUSSION AND CONCLUSION 

We have studied in this work relativistic stars in scalar tensor and f(R) theories that use the chameleon mechanism. 
The behaviour of the scalar field is extremely similar in the two types of models. The main properties are the following. 

Deep inside the star, the scalar field follows very closely the minimum of its effective potential (if it exists). As a 
consequence, if the minimum changes as a function of the radius, the scalar field will closely follow it. If p — 3P < 0, 
which can occur for instance in the central region of very compact stars with constant energy density, there is no 
minimum for the effective potential. It is however possible to find numerical solutions for constant energy stars with 
both w < 1/3 and w > 1/3 regions, but only up to some maximum gravitational potential. We have given qualitative 
arguments to show that one expects tachyonic instabilities if most of the star is characterized by w > 1/3. 

As already emphasized in our previous work Q , our results invalidate the claim that highly relativistic stars cannot 
be constructed in f(R) theories. Numerically, the construction of a relativistic stars in f(R) gravity is a challenging 
task because the scalar field value is extremely close to the singularity in the center of the star. We have also 
constructed relativistic stars in the so-called "cured" f(R) models, which have been advocated to solve a cosmological 
singularity problem. 

We have also illustrated how the screening effect of the chameleon mechanism manifests itself for relativistic stars, 
by computing numerically the effective coupling of the star as well as the post-Newtonian parameter j. 



Note: while this paper was being completed, another paper |37| on relativistic stars in f(R) theories appeared on 
the arXiv, although in a different context. 
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Appendix A 
1. General formulas for f(R) models 
For f(R) models, the potential V, introduced in (|29p . can be written as 



V = Mi 



Rf,R f 
2f 2 ~ 

J ,R 



= M 2 p RqV{x), V(x) 



_ xT'{x)-J : {x) 



2F{x) 2 

where f(R) = RqJ-(x) and x = R/Rq. The associated scalar field is defined by 

<t> = ^|m p lnf k = sJ*-M P \nF'{x) 



(Al) 



(A2) 



Combining (|A1[) and (| A2|) , one finds for the first and second derivatives of the potential V in terms of <f> the following 
expressions: 



dV [2 V-Rf,R 



,R 



d 2 V 



1 



3/ 



RR 



RjM *ff. 

f.R I" 



RR 



(A3) 
(A4) 



2. Starobinsky's models 

Let us now restrict ourselves to the specific class of models introduced by StarobinskyQ, 

R 



f{R) = Ro 



-A 1-1 



Ro 



Substituting this expression into (|A3|) and solving V' = yields the minimum of V corresponding to the asymptotic 
de Sitter solution i?^ = XoqRq. It is convenient to express the parameter A in terms of x^: 



A = 



2[(l + x 2 oc ) n + 1 -1- {n+l)s 



In terms of Xoo, the minimum of the potential is given by 

1 



for 



Voo = V^oo) = -Mf,R x c 



M P In 



(1 + a&Jn+i - 1 - (2n + l)x 2 



{l + x 2 QO ) n + 1 -1- (n + l)a 



At the minimum, the effective squared mass is 



d 2 V 



414 



(a* + 1)" +2 - (2n 2 + in + 1) xt - (n + 2)x 2 c 



M 2 



3na&, [(2n + 1)^ - 1] 



The parameter Xqq must be chosen so that m 2 > 0. 
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Let us now consider the regime where the curvature is much higher than the cosmological curvature Rq, i.e. i>1. 
In this regime, 

f(R) ~R (x-\ + \x- 2n ) , (A5a) 
f'{R) ~l-2\?ix- 2n -\ (A5b) 

f"(R) ~ -J-An(2n + l)^ 2 ™- 2 , (A5c) 
Rq 

and the scalar field behaves like 

Substituting in (|A3[) the above expressions, one finds that the derivative of the effective potential, whose general 
expression can be found in (|13p . is given by 

The minimum of the effective potential is thus determined by 

p- 3P 
Xmin ~ Mp^' 

and exists only if the matter term p — 3P is positive. Substituting in (|A4j) the expressions (|A5|) . one also finds that 
the effective mass at the minimum behaves like 

2 _ -^0 2n+2 

cff ~ 6An(2n + 1) min ' 

In the particular case n = 1, one can invert the relation between <p and R and obtain an explicit formula expressing 
x as a function of <f), which we do not write down because it is rather ugly. Substituting this relation 

A (x 2 - 1) (x a + x) 2 , 

V(x) = i ^ 77 A6 

2(a; 4 + 2a; 2 -2Aa ; + l) 2 

thus gives an explicit expression for the potential as a function of (f>. 
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